% 加载待分类的数据,210*7
X = load('seeds.mat');
X = X.X;
% 分类的个数
K = 3;
% L(θ,θj)变化小于e就终止迭代
e = 10^-4;

% gamrnd 产生伽马分布的随机数,gamma:210*2,表示每个个体属于每一类的概率,相当于γik
gamma = gamrnd(ones(size(X,1), K), 1);
% B = repmat(A, m, n) 将矩阵A复制m*n块，即B由m*n块A平铺而成
% sum(gamma, 2)将gamma按行相加，列数变为1，将gamma每列做了归一化,表示属于第一类和第二类的概率和为1
gamma = gamma ./ repmat(sum(gamma, 2),1 , K);

%  max(A,[],2)
%  是包含每一行的最大值的列向量,m对应每个个体最大的哪一类的概率,labels对应索引,label为1或者2,取每个个体概率大的那一类的索引。
[m labels] = max(gamma, [], 2);
% 在没有经过EM下的原始分类效果图
Clusters =  plot_graph(K, labels, X);

% 原始分类情况下mu, sigma, lambda的初始化
ll = -inf;
% 最大化L(θ,θj)情况下的mu, sigma, lambda的初始化
[mu, sigma, lambda] = M_step(X, gamma);
% 计算L(θ,θj)
nll = log_likelihood(X, mu, sigma, lambda);
disp(['the log likelihood = ' num2str(nll);]);

% 终止条件：L(θ,θj)变化小于e就终止迭代
while (abs(ll - nll) >= e)
    ll = nll;
    % 每迭代一次就更新一下每个个体属于不同类别的概率
    gamma = E_step(X, lambda, mu, sigma);
    % 最大化L(θ,θj)情况下的mu, sigma, lambda的更新
    [mu, sigma, lambda] = M_step(X, gamma);
    % 计算L(θ,θj)
    nll = log_likelihood(X, mu, sigma, lambda);
    disp(['the log likelihood = ' num2str(nll);]);
    % 每迭代一次就更新一下每个个体属于各自的类别，并展示
    [m labels] = max(gamma, [], 2);
    Clusters = plot_graph(K, labels, X);
end








